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; ABSTRACT 

■ We present new analytical distribution functions for anisotropic spherical galaxies. 
I They have the density profiles of the 7-niodels, which allow a wide range of central 

density slopes, and are widely used to fit elliptical galaxies and the bulges of spiral 
galaxies. Most of our models belong to two two-parameter families. One of these pa- 
^ ■ rameters is the slope 7 of the central density cusp. The other allows a wide range 

I of varying radial and tangential anisotropics, at either small or large radii. We give 

■ analytical formulas for their distribution functions, velocity dispersions, and the man- 
ner in which energy and transverse velocity arc distributed between orbits. We also 
give some of their observable properties, including line-of-sight velocity profiles which 

\^ I have been computed numerically. Our models can be used to provide a useful tool for 

i creating initial conditions for N-body and Monte Carlo simulations. 

i-C ' Key words: galaxies: kinematics and dynamics - galaxies: structure. 

o 

to 1 INTRODUCTION 



A galaxy is fully described by a distribution function J-{r, v) wliich specifies the positions and velocities of all of its constituent 
stars and gas. Eddington (1916) showed how to determine the distribution function for a known spherical mass provided that 
the galaxy is dynamically isotropic, that is with the binding energy E being the only isolating integral. Spherical form is of 



5h I course an ideali sation, and at best an approxi mation to the true shape. Isotropy is generally also an approximation to the 
true dynamics ( lUingworth 19771 : iBinnev 1981 ). Our anisotropic models allow us to explore a range of possible dynamical 
behaviour while still retaining the other simplifying approximation of sphericity. 

The step from finding isotropic distribution functions f or spherical galaxies to that of findin g anisotropic distribution 



functions for them has turned out to be surprisingly steep ( Deionghe 19861 : Hunter fc Qian 1993h . Anisotropic distribution 



functions depend also on the modulus of the angular momentum L. Because of the limited supply of analytical models, many 
workers have followed Hernquist (1993) in finding approximate distribution functions by first solving the Jeans' equations for 
the velocity dispersions and then using Gaussians to provide local velocity distributions. Recently Kazantzidis, Magorrian 
& Moore (2004) have shown that this procedure can be hazardous when applied to generate initial conditions for galaxies 
that are strongly non-Gaussian. In such cases, the Jeans plus Gaussian approximation can lead to an initial state that is far 
from equilibrium, so that much of its subsequent development is an artifact of that start. As the importance of numerical 
simulations increases, the issue of controlling numerical errors in them becomes equally important. One can avoid invalid 
initial conditions for N-body and Monte-Carlo simulations simply by using a three-dimensional Monte-Carlo simulator on an 



exact distribution function to derive the initial conditions ( Buvle et al. 200a ) 



Some of the currentl y known anisotropic models have constant anisotropies llvan der Marel et al. 2000l:[wilkinson et al. 20041 : 



Treu fc Koopmans 2004 *). Others are of the Osipkov-Merritt type ( Osipkov 19791 : Merritt 19851 : Hernauist 1990h and are lim 



ited to depending solely on an isolating integral which is a linear combination of E and . Models of this type are isotropic 
at small radii but dominated by radial orbits at large radii. Cuddeford (1991) showed how to modify the Osipkov-Merritt 
algorithm to allow for anisotropy at small radii and a central cusp, while Ciotti fc Pellegrini (1992) constructed composite 
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Osipko v-Merritt mode l s. Models with a greater variety of anisotr opic behaviour are known for special mo dels; the Plummer 
model ( Plummer 1911 : Deionghe 1987 ) and the Hernquist model ( Hernguist 199C1 : Baes fc Deionghe 2002 ). A greater variety 
of analytical systems with different and varying kinematics is clearly desirable. Some have recently been provided by An & 
Evans (2006a). Theirs are for two families of densities; the generalised isochrone family which includes the Hernquist and 
isochrone models as special cases, and a generalis ed Plummer family which includes t he Plummer sphere. Ours differ because 



they are for the densities given by the 7- models (IDehnen 1993 Tremaine et al. 1994 ]_ 



The 7-models include the models of Hernquist ( Hernquist 1990l), Dehnen dPehnen 19931) and Jaffe ( Jaffe 198^ ) as special 
cases. Isotropic distribution functions are known for all ^-models l|Tremaine et al. 199'4 : lBaes et al. 2005 ). The two new families 
of analytical anisotropic spherical models, which we present in this paper, introduce an additional parameter q which allows 
th e anisotropy of the system to be varied. One family is isotropic at large radii, and q/2 gives Binney's anisotropy parameter 
P (jBinnev fc Tremaine 19871 ) at its center. The second family is isotropic at small radii, and —q/2 gives the value of /3 at large 
radii. A greater variety of behaviour can be obtained by combining the two families, though we do not pursue that possibility 
here. 

In Section[2]we review the basic equations needed for our study, and introduce the 7-models and their isotropic distribution 
functions. Before introducing our first new family in Section [S] we first give some surprisingly simple constant anisotropy 
/3 = 1/2 models of which only the 7 = 1 case was known previously. We then show how to construct the first family of 
variable anisotropy models. We give analytical expressions for their distribution functions, velocity dispersions, and energy 
and transverse velocity distributions. Section [3] concludes with a numerical investigation of some line-of-sight velocity profiles. 
Section [4] gives the same properties for our second family of models with its different anisotropic behaviour. The analysis here 
is compact because it makes extensive use of that done in Section [3] for the first family. Section |4] concludes with another set 
of constant anisotropy models; ones composed entirely of radial orbits. They exist for all 7-models with 7^2. We summarize 
our conclusions in Section [S] The Appendices give mathematical details of our derivations, and collect formulas. 



2 BASIC PROPERTIES 
2.1 General formulae 

The isotropic motion of stars and gas in a spherical galaxy is described by a mass distribution function J-[E) which is solely 
a function of the binding energy E. The mass density p is obtained from the distribution function by the integration 



p{^) = A7T T{E)^2{i,-E)dE (1) 



where ip is the gravitational potential. Eddington (1916) solved the integral equation ((T} to get 



diP ^2{E-iP) 



where De denotes differentiation with respect to E. Eddington's solution makes use of the mass density p(V') expressed as a 
function of the potential. Further kinematic information about the system can be obtained once the density is expressed in 
this form. The velocity dispersions for example are given by 



c4 

a^W = a^W = ali-^) = a^o{^) Pii'')di'' ■ (3) 

The distribution function of an anisotropic model of a spherical galaxy function depends also on the modulus of the 
angular momentum vector 



L = r^/v^+vj = rvt. (4) 
Then a double integration is needed to derive the mass density from the distribution function J-{E, L) as 



p{M = 2^ dE ;^'''\> dvl (5) 

Jo Jo y/2{i, -E)- 

It gives what Dejonghe (1986,1987) has named the augmented mass density p('i/),r). Unlike the p(;ip) of the isotropic case in 
equation ([l]), this augmented mass density is not determined by the spatial dependence of p{r) and i/'(r). Instead different 
feasible distribution functions J-{E, L) give rise to different augmented densities for the same mass density. The inversion 
of equation ^ to determin e the distributio n function which corresponds to a known augmented mass density is generally 
unstable if done numerically I Deionghe 19861 ). Velocity dispersions can be derived directly from the augmented density, without 
determining the corresponding distribution function, as 



Completely analytical families of anisotropic ^-models 3 

1 1 

al{tl},r) = crg{ip,r) = -at{ip,r) = — — -/ D,.2[r^ p(tp' ,r)]dili' . (7) 

However, until the distribution function has been determined, one must beware that these dispersions might correspond to 
an unphysical model. 



2.2 Isotropic 7- models 

The 7-models were independently introduced by Dehnen (1993) and Tremaine et al. (1994) and are a Renera lisation of a 
series of spherica l models su ch as the Hernquist model fife rnauist 1990l . 7 = 1), the Dehnen model ( Dehnen 19931 . 7 = §) and 
the Jaffe model i Jaffe ig's^ . j — 2) that are famous for their analytical description of the observed cuspy slopes in elliptical 
galaxies. The dimensionless mass density of these models is given by 



p(r) 



47r r-f {I + r)*-'i '■ 



(8) 



where 7 can have a value between and 3 and determines the growth of the density at small radii while density decays as 
at large radii. The potential of a 7-model is given by 



t/)(r) 



■ 7 



1 



1 + r 



2-7 



for 7 7^ 2. The combination of equations (|8} and ([9]) gives the density as the following function of the potential: 

3_^, — r , -,4 



pW = ^ [1 - (2 - 7)^]^ {1 - [1 - (2 - 7)V]^} 



(9) 



(10) 



Tremaine et al. (1994) used equations ([2]| and (|10p to derive isotropic distribution functions for all 7-models. Baes et al. (2005) 
showed that they can all be expressed in terms of hypergeometric functions as 



(7 - 4) 2F1 1 



2(7 - 1) 2F1 



2-7' 2' 

I ^ 3 

'2-7' 2' 



-;(2-7)i5 +2(7-3) 2F1 1, 



1-73 



o;(2-7)s 



(2-7)£; +7 2F1 1, 



-7 3 



2 - 7' 2 



-7'2 
(2 - 'y)E 



(11) 



The distribution function for the 7 = 2 case of the Ja ffe (1983) model can be reco vered by taking the 7-^2 limit. Then 
the hypergeometric functions become confluent ones (jAbramowitz fc Stegun 19651 ) according to limT,_,2 2^^1(1, a(7)/(2 — 
7); 3/2; (2 — ^)E) = M(l, 3/2, a{2)E). The result is equivalent to that given by Jaffe in terms of Dawson's integrals. 



3 ANISOTROPIC 7-MODELS 

The first of our two-parameter families of anisotropic models is introduced in ^3.21 A preliminary section i )3.1l gives a set 
of constant and radially biased anisotropic models whose distribution functions are much simpler than the isotropic ones 
(|ll|l . These models exist for all 7 ^ 1, and resemble components of the family which is the main topic of this section. The 
distribution functions of this family are given in i^3.2l as infinite series. Convergence requirements generally restrict this family 
to the range < 7 < 2, though particular models for which the series can be summed explicitly, such as those of ii3.2.1l 
exist for larger 7 values. Later sections give the velocity dispersions, distributions of energy and transverse velocity, and line 
profiles of this family of anisotropic models. 



3.1 Simple models for 7 > 1 

An & Evans (2006a) give rules for deriving anisotropic distribution functions which depend on angular momentum via a power 
law. The case of an inverse first power is particularly simple. We first extract an r^^ power from the augmented density (|10p . 
and then apply their algorithm. The result is the one-parameter family of distribution functions 

^(i5,L,7) = |^|l-[l-(2-7)i5]^|'|4-7+ '^—^ -\ ■ (12) 

1^ J I [1- (2-7)£]— J 
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This family exists only for 7^1 because of the second term in the last brace. That term, which becomes large and dominant 
as (2 — tends to its central value of 1, is negative when 7 < 1. The restriction to 7 1 is a simple instance of An & 
Evans's (2006b) cusp slope-central anisotropy theorem because Binney's anisotropy parameter (3 has the constant value of 
1/2 for the models (|12|l . Two examples are 

^(^'^'^) = l^' ^(^'^'2)- 8.3^ (13) 

The first was given by Baes & Dejonghe (2002), while the second can either be obtained directly, or from the 7^2 limit of 
equation (I12|) . 



3.2 A two-parameter family of anisotropic models 

We now construct a family of anisotropic models for different 7-models by writing the mass density ((8)l in the separable 
augmented form of r-~'(l + r)' times a function of ?/' as follows: 



= I, [l-(2-7)^]^{l-[l-(2^7)V'l^} 



3-7 (1 + r) 



47r ^-^ \ k 

fc=0 



2-T 



- ^^E<-)'p-)*i'i:<-')'(:)(T)' '»> 

p=4 fc=0 \ / \ / 

The parameter q allows a range of models for each mass density ((Sjl. The binomial expansion in powers of (2 — 7)^ is valid 
provided that 2 > 7 > because then (2 — 7)1/) G [0, 1]. It begins with p = 4 because the augmented density (I14|l varies as 
ip'^ for small ip, i.e. at large distances. Series converge rapidly except in the central regions where the distribution function 
becomes infinite because of the cuspiness of the mass density there. 

We define Fg{E, L) to be the component of the distribution function which corresponds to the component r~''(l + rYip^ 
of augmented density. It is given by formulas (|B1|I and HB3|) in Appendix [B] The full distribution function is then found by 
summing 

00 

HE,L,^,q) = ^^C(7,P,g)J',''(i?,L), (15) 

p— 4 

where the coefficients C(7,p, g) are defined as 

C(7,P,g) = (-l)''(2 - E(-l)' (fc) (^) ■ (16) 

The distribution function components F^{E,L) do not depend on the parameter 7. This parameter influences the rapidity 
with which the series (|15p converges through the coefficients (7(7, p,g), and especially their dependence on the power (2 — 7)''. 
Hence the closer 7 is to 2, the more rapidly does the series (|15p converge. 

The elementary distribution functions are simpler in certain special cases which are featured in the figures. For q G 
[1, 0, —1, —2], we have 



(27rS)3/2 



r(p) q^ V / 1' g = o,i, r 



(17) 



As we see from equation (|24[) below, q/2 gives the central value of Binney's anisotropy parameter, and hence An & Evans's 
(2006b) cusp slope-central anisotropy theorem restricts q to the range g ^ 7. It is also necessary that q <2 because grows 
as I/"' as 1/ ^ when g > (See equation l|B3[) l. and the integration over in equation ([5} diverges if g J5 2. 
The expansion of the augmented density in powers of tp is simple for the 7=1 Hernquist models, for which 

C{l,p,q) = {~-ir(l^^X p^4. (18) 



p — 4 

We display distribution functions for these models for different g values in Fig. [T] They are plotted there in peri- and 
apocentre space. Circular orbits lie along the lower diagonal boundary, while orbits become increasingly radial as the left 
vertical boundary is approached. Hence densities are larger near there for the radially anisotropic g = 1 case relative to the 
isotropic g = case, but relatively smaller for the tangentially anisotropic g = — 1 and g — —2 cases. 



q=l 
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1 2 3 

pericentre 
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Figure 1. Turning-point diagram for several Hernquist models of the first family of i|3.2l The colour shows how the density of orbits 
varies with pericentre and apocentre. Density declines from a high of white at the very centre, and then through the spectrum from red 
to a low of blue. Note the high red densities of elongated orbits along the left boundary for the most radial q = 1 model, and their steady 
decline with increasing tangentiality as q decreases. 



3.2.1 Explicit q = 1 distribution functions for 7 J? 1 

Equations (|15p and (|17[1 show that the distribution functions for g = 1 have the simple form 
:F{E,L,-y,q = l) = fo{E) + i^. 



(19) 



The two terms correspond to the two components po{ip) and pi(;ip)/r of the augmented density (|14p . That allows us to avoid 
the expansion in powers of xp, and to calculate the two components directly combining the methods of Baes et al. (2005) and 
An & Evans (2006a) respectively as in Sections 12. 21 and 13. II The results, which are similar to but not the same as the earlier 
(fTTI) and ini, are 



ME) 



and 



h{E) = 



87r» 



'2E 



-3(7 - 5) 2F1 ( 1, I (2 - 7)S j + 8(7 - 4) ^ 1, 2 



7 2 



(2-7)S 



2-7' 2' 



-6(7 - 3) 2F1 ( 1, I- (2 - 7)^5) + (7 - 1) 2F1 ( 1, 1-4; I- (2 - 7)S 



2-7' 2' 



|l-[l-(2-7)i5]^}'i5-7 



7-1 



[l-(2-7)iJ] = 



(20) 



(21) 
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Figure 2. Radial velocity dispersion profiles for models of the first family for three 7 values. Dispersions increase with q, and the top 
curve corresponds to the highest allowed value 5 = 7. 



This solution is also restricted to 7 ^ 1, like that of S ection r3.il because otherwise JF becomes negative in the physical range. 
However 7 may exceed 2, because the limitation imposed by the series expansion in used to obtain the solutions in Section 
3.2l no longer applies. The distribution functions are particularly simple for the 7 = 1 Hernquist case for which 

8E^V2E E'* 



T{E,L,1,1) = 



(22) 



3.3 Second-order moments 

An analytical expression for the velocity dispersions of all models can be derived by means of equations (|6]) and ((7]). Since the 
augmented density is separable in r and ip and depends on r as (1 + r)''r~'', the tangential velocity dispersion can easily be 
derived from equation ((Tjl as 

aUr,q)^{l-l^)a'r(r,q). (23) 
Hence Binney's anisotropy parameter is the monotonic function 

for all 7. It has the same sign as q. Systems are more radial than isotropic when q > 0, and more tangential than isotropic 
when g < as seen in Fig. [3l though they all tend to isotropy at large r where /3 — > 0. The condition a'g{r,q) restricts 
q ^ 2. There is no lower limit on q, and the limit q — > —00 gives a system with all orbits circular. 

The simple separable form of the augmented density in equation p4p . combined with equation (O , leads to a compact ex- 

plicit expression for the radial velocity dispersion for all 7 and q in terms of an incomplete Beta function l Abramowitz fc Stegun igGSl ) 
as 



CFr{r,q) 



S^(5,2-27 + g), 



(25) 



cf i Baes fc Deionghe 2OO2I ). Fig. [2] shows the radial dispersion ar{r,q) tending to a common form at large r where all models 
tend to isotropy. Equation (pS)) shows that common form to be l/\/5r. The figures also show that ar{r,q) increases with 
increasing q near the centre as the orbits there become more strongly radial. 



3.4 Energy distribution 

The anisotropic models difler in way in which energy is distributed among their orbits. We label the part of the mass density 
that is contributed by the energy range [E, E + dE] as the energy density J^e- It is given by the inner integral of equation ([SJ 
as 

Jo ^2{xb ~E)- vf 

We define as ^{tp,r,E) the components of energy density which are obtained with the elementary distribution functions 
Fq{E,L) on the right hand side of this equation (|26p . They are evaluated in Appendix]^ Then the energy density for the 
full model is given by the sum 

00 

J'Eitb,r,E) = ^^C(7,P,g)F|,,(V',r-,iJ), (27) 

p— 4 
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Figure 3. The variation of the anisotropy parameter with radius r for the models of i|3.2l It is independent of the model parameter 7. 
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Figure 4. Energy densities for the 7 = 1/2 and 7 = 3/2 models of the first family at spatial radii r = 0.1, r = 0.8, and r = 5. The same 
colour-codings of the values of q are used through each row. The energy E is scaled with its maximum value V(^) that radius. Each 
curve in a panel encloses the same area. That area is the local mass density divided by -0, and varies from panel to panel with both r 
and 7. The singular growth of the g = 7 = 3/2 curve for r = 5 does not occur until E/tp > 0.999. 



with the coefficients C['y,p, q) that were introduced in equation (|16p . The energy density for the simple 7 = g = 1 distribution 
function (1221) is 



64 I E_ 



5/2 



1 - - + — I — 

ijj 7rr \ -i/i 



1/) 



(28) 
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Figure 5. The transverse velocity densities for the 7 = 1/2 and 7 = 3/2 models of the first family at spatial radii r = 0.1, r = 0.8, and 
r = 5. The same colour-codings of the values of q are used through each row. The transverse velocity is scaled with its maximum value 
\/'i'i> at that radius. All the curves in a single panel again enclose the same area; it is now the local mass density divided by y'ixj). 



Fig. |4]shows the energy densities for two values of 7 at three radii: one small, one intermediate, and one large. We chose 
the 7 values of 1/2 and 3/2 to represent the range < 7 < 2 for which our solutions apply. Te is more concentrated at higher 
energy at small r for the more tangential models (smaller q) and less concentrated for the more radial models (larger q). The 
situation is reversed at large r. The reason is that the more radial orbits are less tightly bound than the more tangential ones 
near the center, while the opposite is the case at large r. Fig.l of Dejonghe (1987) for anisotropic Plummer models shows 
the same phenomenon, though our are more concentrated towards high energies than his because our 7-models are more 
centrally concentrated than his Plummer model. Also, unlike his, the differences between our models diminish at large r where 
they become isotropic. 

Whereas Te drops to zero at the upper limit E/tp — 1 for all of the upper row of 7 = 1/2 models, and also for the q ^ 1, 
7 = 3/2 models of the lower row, it becomes infinite as E/ip 1 for the most radial and most cuspy g = 7 = 3/2 model, 
and has a finite limit for the q = 1 model, as in the 7 = 1 case in equation (I28|) . This behaviour is a consequence of the 
{1 — E/^p)^^~''^^^ dependence of ^ for E/ip > r^/(l + r^) given by equation (|B7|I . That dependence is a consequence of the 
L"' dependence of the distribution function for < 2E as given by equation (|B3|l . 



3.5 Distribution of the transverse motions 

The anisotropic models also differ in way in which transverse velocity is distributed among their orbits. We label the part of 
the mass density which is contributed by the transverse velocity range [vt,vt + dvt] as the transverse velocity density J^vf It 
is given by the inner integral of equation ([Sjl , after the order of its integrations has been changed, as 

T.,(i^,r,vt)=4nvt r —£^l^=dE. (29) 
Jo ^2{tp -E)- 

We evaluate the transverse velocity densities F^^^q{^,r,vt) corresponding to the elementary distribution functions Fg{E,L) 
in Appendix!^ The transverse velocity density for the full model is 
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^vt (ip, r, vt) = 



3-7 
4-K 



Vt), 



(30) 



p=4 



with the same coefBcients C{'y,p,q) as before. The transverse velocity density for the simple 7 = 5=1 distribution function 
(I22I) is 



^vtiip,r, Vt) 



7/2 



4-K 



t{i-ty + 



32 

SSvrr 



2^7/2 



^ t 



Vt 



< 1. 



(31) 



Fig. [5] shows the transverse velocity densities for the same set of models as those whose energy densities were shown in 
Fig. m Now J^y^ is generally more concentrated toward low velocities at small r for the more tangential models (smaller q) and 
less so for the more radial models (larger q). The situation is reversed at large r, though differences there are small because 
our models tend to isotropy there. Fig. 4 of Dejonghe (1987) for anisotropic Plummer models shows similar behaviour at small 
radii. The generally higher transverse velocities at small r for the more radial models is due to orbits which are close to their 
turning points and so have transverse velocities which are well in excess of the local circular velocity. 

The L~' dependence of the distribution function for < 2E again has an important effect. It is now seen as Vt 0. 
Equation (|B11[) shows that F^^^q varies as vl~'' for vt/V2^ < l/Vl + r'-^- It tends to as «i ^ for g < 1, has a finite limit 



there for q — 1, and becomes infinite for g > 1. The infinite growth is more prominent for the 7 ■ 
for those of Fig. |4] because of its lower power (-1/2 rather than -1/4). 



3/2 models of Fig. [5] than 



3.6 Observable properties 



The projected density of a 7-model is found from the relation 
rp{r)dr 



E{R) 



\/r2 - 7?2 ' 



(32) 



where we align the z-axis with the line-of-sight, and R — \/ x'^ + is the usual plane polar coordinate, now in the plane of 



the sky. The projected velocity dispersion o-p{R) is given by (jBinnev fc Mamon 19821 : iBinnev fc Tremaine 19871 ) 



aliR) = 



E(i?) 



p2 



r p{r)a'^{r)dr 



7?2 



(33) 



The integrations needed for both (|32p and ((33|) must be done numerically for 7-models. Projected densities of some 7-models 
are shown in Fig. |S] the dependence of the central density slope on the parameter 7 is clearly visible in projection. Projected 
velocity dispersions are shown in Fig. [T] They decrease monotonically with the projected radius R in radial systems, but can 
have a central peak for tangential systems when the tangential component of the velocity dispersion becomes increasingly 
important with increasing R. 

The normalised line-of-sight velocity profile l{vz,R) describes the distribution of along the line-of-sight. It is obtained 
by integrating the distribution function over the velocities vr and in the plane of the sky, followed by a spatial integration 
along the line-of-sight. It is given by 

T.{R)l{v,,R)=2 I f I dvRdv^T{E,L). (34) 



Vr2 - _R2 



It falls to zero at the extremities Vz = ±yj2il)(R). We express the arguments of the distribution function T as 
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log(R) log(R) log(R) 



Figure 7. The projected velocity dispersion profiles for the same set of models as in Fig. |2] 





Figure 8. The normalised line-of-sight velocity profiles at projected radii R = 0.1 and R = 5 for models with three different values of 
7. The same colour-codings are used in each column, i.e. for each 7. 



E = i)-^{vi + vl + vl), L=y^r^v^ + {zvR-Rv,y, (35) 

and integrate using polar coordinates in {vr, v,p) velocity space. The area under the normalised velocity profile is unity because 
integrating JT over all three velocity component gives the density p{r). Hence integrating the right hand side of equation 
l|34|) over Vz gives the projected density E, and hence integrating both sides of equation (|34|) over leads to the result 
/ l{vz,R)dvz = 1. 

Fig. |8] shows normalised line-of-sight velocity profiles for models of the first family for three values of 7, all obtained by 
evaluating equation (|34|l numerically. The most striking feature of these profiles is how much more they vary at small radii 
with the central density slope 7 than with the orbital composition. The profiles of the more tangential models are a little 
narrower at small radii for all 7, and slightly broader at large radii where they are almost isotropic. Other features that stand 
out are the discontinuous slopes at = of q = 1 profiles, and the sharp cusps there of g = 1.5 profiles. They are further 
visible effects of the L"'' dependence of the distribution functions for < 2E. 
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1 

0.5 
P 
-0.5 

-1 

-2-1 12 

log(r) 

Figure 9. The anisotropy parameter /3 as a function of radius r for the second family of models of ^4.11 It is again independent of the 
model parameter 7. 

4 MORE ANISOTROPIC MODELS 

The models introduced in Section [3.21 are all isotropic at large radii. If isotropy is a consequence of mixing, then it is more 
likely that galaxies are isotropic in their well-mixed centres than in their outer regions. We now construct a second family of 
models which are isotropic at their centres but anisotropic at larger radii. We do this by modifying the method of the previous 
section, in a way which allows us to use similar mathematics. We derive the distribution functions of these models in ^4.11 
We give their velocity dispersions, energy and transverse velocity distributions, and projected line profiles in the next two 
sections. Then in a brief final section 8^4.41 we construct some extreme 7-models for which all orbits are radial. 




4.1 Another two-parameter family of anisotropic models 

We obtain these by writing the augmented density in the form of (1 -f r)' times a function of ?/) as 

= _ (2 _ 7)V.]5^ |l - [1 - (2 - 7)^]^}'^' • (36) 

As a result, Binney's anisotropy parameter is now 

2(iT7)' 

and has the opposite sign to q. It is zero at the isotropic centre, and tends to — g/2 at large distances (see Fig. [9]). Systems are 
more strongly radial in the outer regions for negative g, and more tangential for positive q. The parameter q is restricted to 
q ^ —2 by the requirement that ag J5 0, but it is no longer restricted by An & Evans's (2006b) cusp slope-central anisotropy 
theorem because the centre is isotropic. Models with the extreme value q — —2 are purely radial at large distances. 

The augmented density (I36|l varies as 1/)''+'' for small ip. Hence its expansion in powers of now ascends in integer steps 
from that initial term, and has the form 
00 

Pi^^ = ^ ( 1 + r )^ ^ C(7, 4 + g + i, g)^'+'^+' , (38) 
j=o 

for suitable coefficients C. The distribution function is 



HE, L, 7, q) = ^ J2 ^(-Y' 4 + 9 + .7. 1)^'^'^' (E, L). (39) 
3=0 



where Fq[E,L) is defined to be the component of the distribution function which corresponds to the component (1 -f r)'i/;'' 
of augmented density. Formulas for it, which again do not depend on the parameter 7, are given in Appendix IB] 

The evaluation of the coefficients C is now more complicated. The method of the previous section in which we expand 
first in powers of [1 — (2 — 7)^'] and then in powers of ip, yields an expansion in integer powers of 1/), and so is valid only when 
(g -f 4) is a positive integer. For those cases 

C(7, 4 + g + g) = (-1)''(2 - lYT^i-^)' (' t ') (4 }j ■ (40) 



12 Fitter Buyle, Chris Hunter, and Herwig Dejonghe 



q=-2 



q=-i 
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Figure 10. Turning-point diagrams, coded as in Fig. [T] but now for Hernquist models of the second family. The predominance of 
near-radial orbits in the q = —2 model, and their decline as q increases can be seen clearly seen. There is no g = plot because it is the 
same as that in Fig. [T] 



As before, the closer 7 is to 2, the more rapidly does the series (|39p converge. 

More generally, one must expand the augmented density (|36p directly in powers of without any intermediate expansion. 
That leads to an infinite series in ip for 1 — [1 — (2 — 7)'i/']^^^^ the need to multiply the (4 -I- g)th power of that infinite 

series with another infinite series. This procedure does simplify when 1/(2 — 7) is a positive integer. The simplest instances 
are the 7=1 Hernquist and the 7 = § Dehnen models for which 



C{l,4 + q + j,q) = l, c(|,4 + g + j,g) =^^(-i)'(j-fc+l)(j-fc + 2) 



4-f g 
k 



(41) 



Although the infinite series expansions (|39} are limited to the range < 7 < 2 by the same convergence conditions as in 
i )3.2l that limitation disappears when the sum in (|38|l is finite. It is finite for the most radial q — —2 case when 1/(7 — 2) is 
a positive integer, as for 7 — 5/2. Equation (|39p then gives exact distribution functions for strongly radial models with steep 
cusps. 

The elementary distribution functions are again simpler for some small integer values of q. Specifically, equation (|B4|l 
gives 



" (27ri5)3/2 \r(p-i; 



2q 



V^r{p - 1) 

for g = 0, 1, 2, and equation (|B2|) gives 



g(g-i) 

2E ' 2r(p-|)2£ 



+ 



(42) 
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EP f V2E\ ^ r(p + l) „ / 1 2E\ , , 

= 7^(2^ y-ir ) r(p 1- ,) -^- - 2 , - . - 1; -xr) - («) 

for g = —1,-2. Exact distribution functions of the forms T(E,L,-^.^q = 1) = fo{E) + Lfi(E) and J^(E, L,^,q = 2) — 
fo{E) + Lfx{E) + L^f2{E) for g = 1 and g = 2 respectively can be found by the methods of An & Evans (2006a). The results, 
which are similar to those of i ]3.2.1l are given in Appendix IB4I 



4-. 1.1 Explicit distribution functions for Hernquist models 



The fact that each C coefficient is 1 for Hernquist models allows the summations needed for the full distribution function 
to be performed explicitly for the simple functions listed above. Those for q — —1 and q — —2 can be done using Euler's 
formula (15.3.1) of Abramowitz & Stegun (1965) for hypergeometric functions which gives 

l--^] = r(p - 1 



2F1 (^-q- ^,-q;p-q 



A. 

L^J r(-g)r(p-i) 7o 



-i(l_t)P-2(^ 



1 + 



2EtY+2 



1/2 J 



dt. 



(44) 



This integral converges for all positive E, regardless of the magnitude of . Using it, the series summation X/j^o^-^ "^^^^ 
l)x^ = 2/(1 — 2;)'' and a change of integration variable to y = Et, we obtain the compact integrals 



1 

2^ Jo (L2+2i/)1/2 
1 ydy 

'tt^ Jo 



dy 



g = -2. 



(45) 



(i-b+m)-^ 

(L2+2a)372(l_B + y)3 ' 

for the two distribution functions. They are evidently non-negative everywhere, and so physically acceptable. The integrals 
(|45[) can be evaluated by elementary methods as 



T(E,L,l,~l) 
and 

HE,L,l,-2) 
where 

X = 2{1 -E)-L^ 



1 

2^3 1 2x 



7i + yjL^ + 2E - 



2{1- E + 2L^)T{E,L,1,~1) 



(l-E) 



+ 



2X 



y/L2 + 2E - 



L 



{l-EY 



+ 2£ + l| 



VL^ + 2E (l-E)^ 



{1- E + 2L^){^/L^ + 2E - L) 



h 



dy 



2 



arctan 



Vx 



„ {L^ + 2yY/^{l~E + y) 1 



~ arctan -== 



( VzJIfiB + v^) { L - V^) 



X > 0, 

x<o. 



Both distribution functions vary smoothly across the curve 2(1 - £■) = where 



T{E, ^2(1-5), 1,-1) 



7r3^ 1 5 



T{E,^2(l-E),l,-2) = 



7r3^ 1 5 



(1-^)5/2 
1 

(1-^)5/2 



1 



Vl - -E - 1 
{l-E) 



7 



(l-£;)7/2 



(46) 



(47) 



(48) 



(49) 



(50) 



We display distribution functions for these models for different g values in Fig llOl They now show increasing radiality 
with decreasing g, and can be compared with those of Fig[T] The g = case, missing from Fig llOl is the same as that shown 
in Figffl 



4.2 Second-order moments 

The velocity dispersions are now related by 

cjl{r,q). (51) 
and the radial velocity dispersions are 

ol{r, q) = r^(l + r)*+^-^Bi/(i+.) (5 + g, 2 - 27). (52) 

Equation (|52|l . when approximated for large r, shows that the radial velocity dispersion (Tr(r, g) ~ 1/ -^(5 + q)r for large r and 
increases as g becomes more negative. The radial velocity dispersion near the centre behaves differently according to whether 
7 is greater or less than 1. It has the form (Jr(r, g) ~ ^ r2-''/[2(7 — 1)], and (Tr(r, g) ~ ^ (— r In r) in the limit 7 = 1, which 
is independent of g when 7 ^ 1. It depends on g for weaker cusps with 7 < 1. Then (Tr(r, g) ~ ^r'^Bi^ + g, 2 — 27), and is 
larger for the more radial models even near the centre as Fig. [11] confirms. 



gr 



2(1 + r) 
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Figure 11. The radial velocity dispersion profiles for three 7 values for the second family, and for q = —2, —1, 0, 1, 2. The highest curve 
is for the lowest value of q and the others follow in sequence. 
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Figure 12. Energy densities, as in Fig.|4l but for now for models of the second family. 



4.3 Other properties 

We define energy and transverse velocity densities ^[•ij^,r,E) and F^_^^q(^,r,vt) corresponding to elementary distribution 
functions Fq{E,L) in the same way as we defined F^ ^{ip,r,E) and Fy^^q{il),r,vt) in ^3.41 and 13.51 The energy and transverse 
velocity densities are then given by sums of the same form as those of equation (|27p and pop but now with barred components. 

Energy densities for models of the second family shown in Fig. [12] show the same relative differences between more 
tangential (now larger q) and more radial (smaller q) as do the models of the first family shown in Fig. U The differences 
are now small at small r because models of the second family are isotropic at their centres. Differences appear as the radius 
increases. The most radial q = —2 model remains strongly peaked even at high E/tp, while the peaks of the other models 
move to lower relative energies with increasing radius and decreasing radiality, as in Dejonghe's Fig. 1. All J^e curves now 
drop to zero at the right limit E/tp = 1 because the distribution functions are not singular as L 0. Our J^e are again more 
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Figure 13. The transverse velocity densities, as in Fig. |5] but for now for models of the second family. 




concentrated towards high energies than those of Dejonghe (1987) because our densities are more centrally concentrated than 
his. 

The transverse velocity densities in Fig ll3l are also similar at small radii. This distribution remains strongly peaked at 
low relative Vt for the most radial q — —2 model at large radii, while the peaks of the other distributions become increasingly 
central with increasing radius and increasing tangentiality. Dejonghe's Fig. 4 shows similar behaviour at large radii. 

Projected velocity dispersions for the second family are shown in Fig. 1141 All now peak at an intermediate projected 
radius R. Unlike the unprojected radial velocity dispersions shown in Fig. 1111 the projected ones become larger near the 
centre and smaller in the outer regions as the radiality increases. 

Fig-Elshows normalised line-of-sight velocity profiles for models of the second family for three values of 7 and at the same 
radii as those of Fig. [8] They too vary more at small radii with the central density slope 7 than with the orbital composition. 
The profiles of the more tangential models are now more noticeably narrower at small radii and broader at large radii than 
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Figure 15. The normalised line-of-sight velocity profiles for models of the second family for the same R and 7 values as Fig. |8] 



the models of the first family. Now that the distribution functions are isotropic at their centres, their velocity profiles peak 
smoothly a.t Vz = 0. 



4.4 Models with all radial orbits 



Models with all radial orbits have distribution functions of the form 



HE,L,'y) = f{E)S{L'). (53) 
With this form, equation (O becomes the following Abel integral equation for f{E) ( Richstone Sz Tremaine 1984h : 



r p = 2n 



f{E)dE 



(54) 



/o ^2(^ - E) 

Physical solutions are possible only if limr-»o p > 0, because (|54|l gives r^p as a weighted integr al of f(E) over the whole phys i- 
cal range of E, and so a zero value would imply that there are unphysical negative values of f{E) (jRichstone fc Tremaine 19841 ). 
The solution of equation (|54p is obtained using the same t-substitution as in Baes et al. (2005) §5.2. It is positive everywhere, 
and hence always physical, because it is an integral with a positive integrand. It can be expressed in hypergeometric functions 
as 



■7. 



4n^ 



'2E 



7 - 2 - 2(7 - 3) 1 



1 



7 - 



2' 2' 



(2 - 7)£; + (7 - 4) 2F1 1 



7 - 



2' 2 



(2-7)i? 



(55) 



This expression provides analytical 7-models for the range 2 ^ 7 < 3, and adds to the only previously-known analytical model 
with all radial orbits (Fridman & Polyachenko 1984, Binney & Tremaine 1987, problem 4-19). The hypergeometric functions 
are elementary for the strong cusp case of 7 = |, as they are for that case in Tremaine et al. (1994), and give 



/2E5{L^ 



+ 



+ 



15 



2(2 + E) 4(2 + £)2 4:{2 + E)s 



21n[yTT| + y| 
VE{2 + E)y2 



15 

4(2-H£)2 



(56) 
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5 CONCLUSIONS 



We have provided a wide range of analytical anisotropic stellar dynamic distribution functions for the widely used 7 density 
profile. The simplest are two one-parameter families with constant anisotropy, those of i]3.1l for 7^1 for which Binney's 
parameter /3 = 0.5, an d the all-radial models of i |4.4l for 7^2 for which P — 1. Only the 7 = 1 member of the first family 
was known previously i Baes fc Deionghe 2002h . 



The two two-parameter families of i l3.2l and g4?T]provide a much wider range. They have a second parameter q in addition 
to 7. It controls their anisotropy. The first family is isotropic at large distances, and q gives the value of 2/3 at its centre. The 
second family is isotropic at its centre, and q gives the value of —2/3 at large distances. A greater variety of behaviour can 
be obtained by combining members of the families, as in Dejonghe (1989) and Baes & Dejonghe (2002), though we have not 
done this. 

Finite expressions are given for some of the distribution functions. These include those in H3. 2. II for models of the first 
family with g = 1, and in Appendix [B4] for tangentially biased models of the second family with q = 1 and q — 2. Other more 
restricted cases are the radially biased Hernquist models of the second family of ii4. 1.11 with q = —1 and q = —2 (most radial), 
and the strongly cusped q = —2, 7 = 2-1- 1/A'^, A'^ an integer, cases identified, though not studied, in ^4.11 More generally, 
distribution functions must be found by summing series of ordinary or generalised hypergeometric functions using formulas 
given in Appendix [B] The coefficients of those series depend on powers of (2 — 7), and the series converge rapidly for 7 close 
to 2 when few terms are needed. Conversely, long series are needed to generate the 7 = 0.5 plots shown in the figures. It is 
not difficult to compute hypergeometric functions; some advice on how to do so is given at the end of Appendix lAl 

The first family of models, which are radially biased at their centres when g > 0, provide an interesting example of 
the cusp slope-central anisotropy theorem of An & Evans (2006b). That theorem states that the central value of f3 can not 
exceed half the density slope 7. In fact our models show that their theorem is somewhat more widely applicable than their 
discussion of it allows. Their proof assumes that the distribution function has a Laurent expansion of a specific form in the 
whole neighbourhood of L = 0. Our equations (|15|l and (|B3[l show that our distribution functions have this form but not in 
the whole neighbourhood of L = 0, only in the region < 2E. A different form, given by equations (|15|l and (|Bip applies 
in < 2E < , which includes a part of the neighbourhood of L = near E = 0. Nevertheless their theorem is still valid, 
and their proof of it requires no essential change, because a Laurent expansion of the form they assume is valid near the cusp 
and its component that is proportional to L^' dominates the behaviour there. The wholly radial models of i^4.4l with jS — 1 
everywhere respect the singular g — > 2 limit of their theorem because they exist only for 7^2. 

Besides distribution functions, we have given and displayed velocity dispersions, energy and transverse velocity densities, 
and observable properties, including line-of-sight velocity profiles. The singular central behaviour of the radially biased models 
of the first family is evident in the g ^ 1 curves in Fig. [l] Fig. [5] and Fig. [S] The latter are particularly significant since 
they plot a directly observable quantity. The 7-models are commonly used to describe the kinematics at the cores of giant 
elliptical galaxies, both to model observations and to derive initial conditions for either N-body or Monte-Carlo simulations. 
Because of the past lack of analytical models, many workers have used Hernquist 's (1993) simplified procedure to generate 
approximate distribution functions. Kazantzidis, Magorrian & Moore (2004) have shown that this procedure can be hazardous 
when applied to galaxies that are strongly non-Caussian, as it can then lead to an initial state that is far from equilibrium. 
We have identified cusped galaxies with /3 ^ 0.5 at their centres and whose line-of-sight velocities are cusped as instances 
of non-Gaussian behaviour. The wide variety of analytic distribution functions of this paper provides alternatives to the use 
of Hernquist 's method. Our distribution functions are in exact equilibrium. They can be coupled with a three-dim ensional 
Monte-Carlo simulator to provide valid initial conditions for N-body and Monte-Carlo simulations i Buvle et al. 20061 ). and so 
avoid the subsequent development being infiuenced by artifacts of the start. 
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APPENDIX A: EVALUATION OF DISTRIBUTION AND RELATED FUNCTIONS 

We defined Fg{E,L) in section [3.21 10 be the distribution function which corresponds t o the elementar y augmented density 



r ''(1 + rYip^ . We find it by taking a Laplace-Mellin transform in E and L — rvt as in jDeionghe 19861 ). and obta 



CE-.cML^^{J'iE,L)} = 



2^/2 



v(3-/3)/2 



(27r)»/2 r(l-/3/2) 



r{p + i) r(/3 - g)r(-/3) 



r(-?) 



(Al) 



Although the Mellin transform of r^''(l + r)' does not exist for non-integer q > 0, the final distribution function which we 
obtain for g < is valid in a continous interval of q, and remains valid for g > by the principle of analytical combination; 
it suffices to check that the proper limits can be taken in Equation (|B3|) in order to check the existence of our results. The 
Laplace transform with respect to a is easily inverted, and we are left with a Mellin inversion integral 

-,fl/2 



F^{E,L) 



r{p + i)EP- 



-3/2 



(27r)5/22'?r(-g) 2ni 



)r(-f + f)r(i 



2^2' 



0Q — ioo 



r(i 



2E 



d 



(A2) 



for Fg{E,L). We have here used applied the duplication formula for the Gamma function ijAbramowitz fc Steeun 19651 ) 
equation 6.1.18 to expand the product r(/3 — q)T{—P). The inversion of the /3-integral in equation HA2p gives, 
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(A3) 



by matching to the definition 5.3.1 in Erdelyi (1953) of the Meijer G function G33. A similar analysis for the elementary 
augmented density (1 4- rY^ji^ gives 
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with a different Meijer G function. Yet other Meijer G function are needed for the energy and transverse velocity densities 
J^E and jr„j . Applying the operator on the right hand side of equation (f26)) to equation (|A2[) and then inverting the ^-integral 
gives 



Applying the operator on the right hand side of equation (|29|) to equation (|A2|) and then inverting the /9-integral gives 



1 

2 

1 _ £ 

2 2 



F^,,,{E,L) = 



r{p + i)vt 

29+i7rr(-(?) 



2\P-1 



(A5) 



(A6) 



The corresponding barred quantities differ from these quantities in that the first two columns of coefficients, which are the 
same in equations (|A3|l . (|A5|l . and (|A6|l . are replaced by those of equation (IA4[) . 

All the Meijer G functions which arise in this paper can be expressed as the sum of two generalised hypergeometric 
functions, using formulas 5.3.5 and 5.3.6 of Erdelyi (1953), as 
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for X < 1. The coefficients (ai, 02, 61, 62) here are (1, i, -f , I - f ) for the models of SjSJand (l, i, -f , i - f ' 
quantities for the models of i]4.1l . Certain coefficient differences are the same for both sets, and those differences have been 
used in equations (|A7|I and (|A8I) . The coefficients (03,63) are (p — fiO) ^^o'' distribution functions, (p — ^'~^) f^'' energy 
densities, and (p, 0) for transverse velocity densities as in equations (|A3p . (|A4[) . (|A5|) . and (|A6[) . Equations (IA7[) and (IA8I) 
simplify when the argument of one of the denominator Gamma functions is zero or a negative integer. Then that Gamma 
function is infinite and one term disappears. 

The 3F2 generalised hypergeometric functions reduce to polynomials if any of the first three coefficients are a negative 
integer, and are 1 if any of those coefficients is zero. They reduce to the simpler 2F1 hypergeometric functions if any of 
the first three coefficients is the same as the fourth or the fifth. Otherwise they can be computed by summing their series 
expansions when their final argument is less than 1 in magnitude, as it is for the two separate formulas HA7|) and (|A8|I . and 
in their applications given below. There is no singularity at intermediate cases of x = 1, because the final arguments of 
the hypergeometric functions are then —1. The programme given in §6.12 of Press et al. (1992) for computing ordinary 2F\ 
hypergeometric functions can easily be extended to the 3F2 generalised case. Series with large a coefficients converge slowly, 
and Press et al.'s idea of combining series summation with integration of a differential equation is then helpful. Generalised 
hypergeometric functions are also available in mathematical software packages such as Maple or Mathematica. 
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B2 The energy distribution 
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B3 The distribution of the transverse motions 
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B4 More exact distribution functions 

The 7-models have distribution functions of the form !F{E,L) = fo{E) + Lfi{E) where 
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where 



Ho) = 2F1 (^1, ^; (2 - j)Ej , t = [1 - {2 - j)E]^ . 

They also have distribution functions of the form J^{E, L) = fo{E) + Lfi{E) + f^iE) where 
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